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Abstract 

We revisit the problem of quark production in high energy heavy ion 
collisions, at leading order in a a in the color glass condensate framework. 

In this first paper, we setup the formalism and express the quark spectrum 
in terms of a basis of solutions of the Dirac equation (the mode functions). 

We determine analytically their initial value in the Fock-Schwinger gauge 
on a proper time surface Q s to C 1 , in a basis that makes manifest the 
boost invariance properties of this problem. We also describe a statistical 
algorithm to perform the sampling of the mode functions. 

1 Introduction 

In heavy ion collisions at ultrarelativistic energies, such as those performed at 
the RHIC or the LHC, the bulk of particle production originates from soft gluons 
that carry a small fraction of the projectile longitudinal momentum [1, 2, 3, 4, 
5, 6]. Because of the infrared singularity present in the emission probability 
of soft massless gluons, the occupation number of these gluons increases as an 
inverse power of the longitudinal momentum fraction x, according to the BFKL 
evolution equation [7, 8]. When it reaches values of order l/a s , non-linear 
processes such as recombinations become important and tame the growth of the 
occupation number - a phenomenon known as gluon saturation [1], 

By virtue of this large occupation number, the dynamics of these soft gluons 
is essentially classical, but non-perturbative because highly non-linear [9, 10]. 
The Color Glass Condensate (CGC) effective theory [11, 12, 13] provides an 
organization principle for this regime of strong interactions, and a calculational 
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framework for computing observables relevant to hadronic or nuclear collisions 
involving such densely occupied projectiles. 

In the study of heavy ion collisions, the CGC has been applied to calcu¬ 
late the gluon yield at early times [14, 15]. At leading order (tree level), this 
amounts to solving the classical Yang-Mills equations with light-cone currents 
representing the fast color charges of the two projectiles. At next-to-leading 
order (one loop) [16], one can extract the terms that contain logarithms of the 
collision energy and show that they can be absorbed into the renormalization 
group evolution -according to the JIMWLK equation [17, 18] - of the probabil¬ 
ity distribution of the above currents. 

The CGC framework can also be used in order to study the production of 
quarks in heavy ion collisions. In this framework, the light-cone color currents 
couple only to gluons (because gluons are the dominant constituents of high 
energy hadrons or nuclei), and quarks are produced indirectly from the gluons 
by the process gluons —► qq. Thus, the quark spectrum is one order higher in a s 
than the gluon spectrum. Equivalently, one may say that the quark spectrum 
is a 1-loop quantity while the gluon spectrum is a tree level quantity. It is 
well known that the single inclusive quark spectrum can be expressed in terms 
of a basis of solutions of the Dirac equation, with a color background field 
that corresponds to the LO gluons (i.e. the classical solution of the Yang-Mills 
equations). The choice of this basis of solutions is not unique. When they 
are chosen in such a way that they coincide with the free spinors v s (k)e lk ' x 
or u s (k)e~ lk ' x in the remote past, they are often called mode functions in the 
literature, and we will also adopt this terminology in the rest of this paper. 

This approach can be used to study the production of fermions or scalar 
particles in any situation where the production is due to a classical background 
field [19, 20, 21, 22, 23]. In the context of heavy ion collisions, this formulation 1 
has been used first in studies of electron production in nuclear collisions [29, 30] 
(although this is a pure QED process, its treatment in the “equivalent photon” 
approximation is very similar to the CGC), and later in a computation of quark 
production in heavy ion collisions [31, 32]. This earlier work was limited in a 
number of ways: (i) the basis of mode functions that was used was expressed 
in terms of proper time and the Cartesian coordinates x, y, z , making the boost 
invariance of the problem highly non-obvious, (ii) the sum over the modes was 
restricted to a subset of all the possible modes, and (iii) the resulting quark 
spectrum may be contaminated by spurious lattice doublers at high momentum. 

The goal of this work is to revisit this study in order to overcome all these 
limitations. In this first paper, we first obtain a new basis for the Dirac mode 
functions, that naturally depend on the proper time r, on the rapidity y and on 
the transverse position x± . These mode functions are indexed by the transverse 
momentum k± and a wave number v which is the Fourier conjugate to 77 , making 
them very convenient for a lattice implementation where the grid covers a fixed 
range in 77 . In order to improve the sampling of the mode functions, we use 

1 In the case of proton-nucleus collisions, or in any situation where it is legitimate to expand 
in powers of the color sources of one of the projectiles, a more direct approach is possible, 
that leads to analytical results at leading order [24, 25, 26, 27, 28]. 
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spinors that are random linear superpositions of all the possible mode functions. 
A proper choice of the distribution of the random weights ensures that the 
exact result is recovered in the limit of infinite statistics. With finite statistics, 
this procedure provides a straightforward way to estimate the statistical errors. 
Numerical results based on a lattice implementation of this framework will be 
presented in a forthcoming paper. 

The contents of the paper is the following : In the section 2, we briefly 
remind the reader of the Color Glass Condensate and of the expression of the 
quark spectrum in this framework. We also show in this section how to choose 
a basis of mode function that makes boost invariance manifest. In the section 
3, we present a statistical method to sample the modes, and derive the formula 
for the corresponding statistical errors. The initial value of the mode functions 
on the forward light-cone (i.e. just after the collision of the two nuclei) is 
derived in the section 4. The section 5 is devoted to concluding remarks. A few 
appendices collect more technical material. The derivation of the expression of 
the quark inclusive spectrum in terms of Dirac mode functions is recalled in 
the appendix A, and an alternate derivation following more closely standard 
Feynman perturbation theory is presented in the appendix B. The appendix 
C discusses a technicality in the derivation of the initial value of the mode 
functions, and the appendix D is devoted to the study of a conserved inner 
product between the mode functions. We make an extensive use of this inner 
product in order to properly normalize the mode functions, and as a consistency 
check at various stages of the calculation. In the section E, we use the QED 
version of the mode functions derived in the section 4 in order to recover the 
electron production amplitude in the collision of two electrical charges. 

2 Quark yield in the CGC framework 

2.1 Color Glass Condensate 

The Color Glass Condensate framework is an effective theory that can be used 
to study the early stages of heavy ion collisions, summarized by the following 
Lagrangian density, 

c = + J£) +MHD - m)i/> , (1) 

written here for one family 2 of quarks of mass m. The color charge content 
of the incoming nuclei is described by the two currents Jf 2 > whose supports 
are restricted to the light-cones, ./{'" oc S(x~) and oc <5(:r + ), in a collision 
at very high energy. These currents fluctuate event-by-event, with a Gaussian 
probability distribution in the McLerran-Venugopalan (MV) model 3 [9, 10] that 

2 As long as we do not include the effect of virtual quark loops, we can consider one quark 
family at a time. 

3 At very high energies, this distribution evolves according to the JIMWLK equation and 
will become non-Gaussian. The MV model may be viewed as a model of initial condition for 
the JIMWLK evolution. 
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we use in this paper. If one is interested in the production of quarks in a given 
collision, one could draw randomly one configuration of Jf 2 ’ an( i n °t perform 
an average over these currents 4 . 

In the gluon saturation regime, the currents J[ 2 are inversely proportional 
to the gauge coupling, 

K 2 ~ - g • ( 2 ) 

For this reason, gluonic observables at leading order are expressible in terms of 
a classical color field that obeys the Yang-Mills equation with the source J 1 + J 2 , 

[D„, F^] = J?+J% . (3) 

One should in principle also impose the covariant conservation of the current 
[Dp, Jf + J%] = 0. This constraint becomes trivial in the Fock-Schwinger gauge, 
x + A~ + x~ ~A+ = 0, since it ensures that the gauge potential vanishes on the 
support of the current. We adopt this gauge in the following. Furthermore, 
for inclusive observables, one can prove that this equation of motion must be 
supplemented by a retarded boundary condition [33], such that the gauge field 
vanishes in the remote past, thereby making this classical solution unique. In 
the saturation regime, this classical gauge field is strong, of order ~ 1/g. 


2.2 Inclusive quark spectrum 

In the CGC framework, the fermions do not couple directly to the currents Jf 2 , 
but only indirectly through the gauge field that appears in the covariant deriva¬ 
tive in the Dirac operator 0 — m. Therefore, the natural order of magnitude of 
the spinors is if) ~ 1, in accordance with the fact that the occupation number of 
fermions is bounded by unity. Thus, observables that contain quark fields are of 
higher order in the gauge coupling. For instance, the g 2 power counting for the 
quark spectrum at LO is the same as that of the gluon spectrum at NLO: both 
are 1-loop quantities, the only difference being the nature of the field running 
in the loop (quark versus gluon). In a fixed background color field, the quark 
spectrum at LO is given by the following formula 5 , whose derivation is recalled 
in the appendix A : 


2w p 


dN q 

d 3 p 


1 y- I" d 3 k 

Wf ^ J (27r)32tc fc ,°4+oo 

’ a,b ’ 


{^b\^ksa) x o 


2 


(4) 


4 Note however that the theoretical basis for doing this is not very robust, and becomes 
inconsistent beyond leading order. Indeed, for fixed loop corrections to observables 

contain unphysical logarithms of the longitudinal momentum cutoff that separate the gluon 
modes that are described by the sources J^ and those that are described as gauge fields. 
These logarithms can be absorbed into a redefinition of the probability distribution W[J] of 
these sources (that must now evolve according to the JIMWLK equation). But this procedure 
only works if one performs an average over the sources. 

5 This formula is true to all orders in the currents J£ and . If one expands it to lowest 
order in these currents (i.e. dN q /d 3 p oc 0((Ji J 2 ) 2 )), one recovers the standard result for the 
process gg —> qq with off-shell incoming gluons, derived in the framework of fc^-factorization 
in refs. [34, 35] (see ref. [24] for this comparison). 
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(u!k = \/k 2 + to 2 ) where ipp^b a f ree positive energy spinor of momentum p, 
spin er and color b (since quarks live in the fundamental representation of the 
gauge group SU(N C ), this color index runs from 1 to N c ). In the absence of 
background field, these spinors are given by 6 , 

= U °(P) e ~ lP ' X (P° = Up) ■ ( 5 ) 

However, it may also happen that the gauge fields at x° —> +oo evolve into 
a nonzero pure gauge configuration. In this case, the above spinor should be 
replaced by a color rotated one: 

(*) = < (p) (*)<W' e~ ip x , (6) 

where SU(N c ) is the SU(N c ) matrix defining the pure gauge background. 

In contrast, ipksa i s a spinor that has evolved over the background color field, 
starting at x° = — oo from a negative energy free spinor of momentum k and 
spin s : 


{U/) x - m) ip ksa (x ) = 0 , lira ip ksa {x) = v s {k)e lkx . (7) 

x u —y —oo 

Note that the subscripts a , b refer to the initial color of the quarks. The 
color they carry at the point x is not written explicitly, and is encoded in 
the N c (color) x 4 (Dirac) components of the spinors. 

The inner product ( • | • ) 0 that appears under the integral in eq. (4) is 
defined by 

{^\x) x o = J d 3 x^(x°,x)x(x°,x) . ( 8 ) 

(In the product ip'Xi a ^ the unwritten color and Dirac indices are contracted.) 
The properties of this inner product are studied in detail in the appendix D. 
In this appendix, we also use its conservation as a consistency check of the 
results that will be derived in the section 4. Note that the formula for the quark 
spectrum requires that one takes the limit of infinite time. As we shall discuss 
later in this section, this is also a requirement for the quark spectrum to be 
gauge invariant. 

Eq. (4) is the expression for the fully inclusive spectrum that we are going 
to use in the rest of this paper. The virtue of this formula is that it reduces 
the calculation of a one-loop graph in a background field to solving a (linear) 
partial differential equation with retarded boundary conditions. Even if this 
can be done analytically only for very simple backgrounds, this problem can in 
principle be tackled numerically for completely general backgrounds. 

6 In this equation and in the rest of this paper, we write explicitly only the indices that 
characterize the initial value of the spinor. A more complete notation would read : 

where a is the Dirac index and b' is the color of the quark at the point x. 
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Note that in eq. (4), the spectrum is summed over all the possible final 
states and over the spin of the tagged quark. In order to obtain the polarized 
spectrum, one simply needs to remove the sum over the spin index a. This 
formula also contains sums over the colors of the initial and final fermion. These 
sums should not be undone, as the spectrum of quarks with a given color has 
no gauge invariant meaning, k and s can be viewed as the momentum and 
spin of the antiquark that must be produced along with the quark to satisfy the 
conservation of the flavor quantum number, as reflected by the initial condition 
for the spinor if>ks in the remote past. 


2.3 Boost invariance 

2.3.1 Change of coordinates 

Since collisions in the high energy limit are invariant under boosts along the 
longitudinal axis', it is convenient to trade the longitudinal components of the 
momenta p z , k z in favor of the corresponding rapidities y v and y^. Likewise, the 
proper time r and spatial rapidity y are more suitable than x° and z to map 
the space-time: 

in(^) • (9) 

Besides the obvious change in the measure dp z /co p = dy, one must alter the 
definition of the inner product so that the integration is on a surface of constant 
r (instead of a constant x°), 


{^\x) T = r J d 2 x ± diyip\T,x ± ,r])e x{T,x ± ,y) . 


( 10 ) 


(See the appendix D.) When doing this, eq. (4) becomes 


dN n 


dy p d 2 p± 87 t(27t) 3 
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0 + 

p ± y p ab \ ^k ± y k saj T 


\i> 


), • (ID 


The boost invariance of the problem implies that the inner product depends 
only on the difference of the rapidities y p — yu- After integration over y^, the 
resulting quark spectrum is independent of the rapidity y p . 


2.3.2 Boost invariant spinors 

The boost invariance can be made manifest at the level of the spinors ipp^ VpCrb 
and ibT themselves. Even when the background field is invariant under 
boosts in the z direction, these spinors depend separately on the momentum 

7 Here, we are disregarding the small-# evolution of the color sources in the incoming nuclei. 
This effect would break the boost invariance of the problem due to gluon loop corrections, 
and make the quark spectrum depend on rapidity on scales Ay ~ a J” . 
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rapidity y and on the spacetime rapidity y. This can be trivially seen on the 
vacuum spinors, whose rapidity dependence can be made explicit as follows 


Kly P °b(T,V,x±) = e^ 7V u a (p±,y = 0) e -^prcoBh( Wp -,) e *x-x > {12 ) 

where M p = p 2 + m 2 is the transverse mass. To turn the prefactor into a 
function of y p — y, it is convenient to define transformed spinors as follows, 

^k ±ysa = y-re^W t/; k±ysa . (13) 


The factor \Jr has been introduced for later convenience. After this transforma¬ 
tion, the new spinors ipk±ysa are boost invariant, in the sense that they depend 
on the spatial rapidity 77 and on the momentum rapidity y only through the 
difference y — y (provided that the background field does not depend on y). 

The boosted spinors introduced in eq. (13) also offer the advantage of obeying 
a simpler form of the Dirac equation where rapidity does not appear explicitly 
in the coefficients. In order to see this, first note that 



Then, multiply this operator on the left by exp( 
gives 


e 


47 7 


7% + 7 3 <9 3 



| 7 ° 7 3 ). A simple calculation 


-An 


o-H 7 


(15) 


From this observation, we conclude that the modified spinors if) obey the fol¬ 
lowing Dirac equation : 


i 




T 


D r 1 + 7 'A 



ip = 0 . 


(16) 


One can see that the coefficients of this equation are independent of the 
rapidity y when the background field is boost invariant (so that there is no 77 
dependence hidden in the covariant derivatives). In terms of the boost invariant 
spinors defined in eq. (13), the inner product on a constant proper time surface 
takes a particularly simple form, 


Wx) T = 


d 2 x±dr] (t , x±,y) %(r, x±,y) . 


(17) 


2.3.3 Mode functions in the v basis 

Another useful transformation is to go from a basis where the spinors have 
a definite momentum rapidity y to a basis where they have a fixed Fourier 
conjugate v to the space-time rapidity y, 

'lpk±ysa ^ Ipkj-i'sa = f dy € ^ ^Pkj^ysa • ( 1 ^) 
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When the background field is boost invariant (i.e. independent of 77), ipk ± ysa 
depends only on y — 77 and the spinor ipk ± vsa in the new basis has a trivial 77 
dependence in exp (ivrj): 

1pk±i'sa('T, Xj_, 77) — G ^ 1pk±vsa (+ *^-L) ■ (19) 


v is a conserved quantum number and the 77 dependence of these spinors is not 
altered by their propagation over the background field. Moreover, the Dirac 
equation obeyed by these spinors is effectively 2 + 1 dimensional, 


i 



D r + i 



T 




gAy) + 7 ! A 



'lftk±Lfsa 0 7 


( 20 ) 


since the 77 dependence can be factored out. 

When calculating the inner product of two such spinors (see eq. (141)), the 
integration over 77 trivially yields a delta function, 


Pp±vcrb sa) T — 2'7T(5(Z' / V ) Id X_\_ Xj_) , (21) 


— ^Pj-i '&b j_ v 1 s a n 


where we denote by [ • | • ] the “reduced” inner product that remains after one 
has factored out the delta function. In this basis, the quark spectrum is given 

by 


dN q 

dy p d 2 p ± 


87t(27t) 3 


E 

7, S = T,4- ‘ 


d 2 k l dv 
(2-7t) 2 27T 


lim 

'—>•+00 


'ijj 0+ 


p ± uab \^ k _l vsa 


( 22 ) 


In this formula, it is tempting to ignore the limit r —> +oo and to interpret 
the resulting expression as the quark spectrum at the finite proper time r. 
One should however consider such a generalization with caution, since it is not 
possible to rigorously define asymptotic states at a finite time. 


2.4 Gauge invariance 

Under a local gauge transformation, a spinor ip(x) is transformed as follows 

ip{x) —> O {x)ip(x) , (23) 

where fl(x) is an SU(N C ) matrix. Since the inner product that enters in the 
quark spectrum given by eqs. (4) or (22) is local, it is gauge invariant provided 
of course that the spinors ip- and 0 O+ 1 are gauge rotated consistently. 

In eq. ( 6 ), we have already indicated that the spinors ■0 O+ 1 should be obtained 
from the free solutions in a null background (5) by an appropriate color rotation, 
if the background field at the time of quark measurement is a nonzero pure 



gauge. The square of the inner product that appears in eq. (4) can be written 
as 


(V'paft ^ksajx 0 


d 3 xd 3 y e ip ' iy - x) 


^ksai^’V) «< x ( p ) 


n(y)rt{x) ul(p)ip ksa (x 0 ,x) ■ (24) 


For the time being, let assume that the background field is a pure gauge at the 
time a : 0 where the quarks are being measured 8 . The factor fI(y)fT(a:) depends 
on this pure gauge, and can be obtained by a Wilson line between the points x 
and y, 

fi(y)f2 t (o:) = C/ 7 (y, x) = P exp (ig j dz^A^z)^ , (25) 

where 7 is a path from x to y in the hyperplane of fixed time a; 0 . Thus, in 
practice one would calculate 


0/^6 I Vw) 


d 3 xd 3 ye zp(v x) ip ksa (x°,y) u a (p) 

X t/ 7 (y, x) uUp)ipksa( x °’ x ) • ( 26 ) 


When the background field is a pure gauge, this does not depend on the 
path 7 chosen between x and y. If one extends the use of eq. (26) to a situation 
where the background field at the time x° is not a pure gauge, one still obtains 
a gauge invariant result, but there is now an ambiguity due to the choice of the 
path. Indeed, Wilson lines U-y and Uy evaluated on two different paths differ 
by a Wilson loop, 

U 7 U^, = P exp (ig j) dz p A^z)^ (27) 

defined over the closed loop qUy ' -1 made of the path 7 followed by the reverse 
of the path 7 '. This Wilson loop measures the chromo-magnetic flux across the 
closed loop, and is therefore equal to the identity only if the background field is 
a pure gauge. One expects this irreducible ambiguity to decrease with the quark 
momentum. On the one hand, in the spectrum of quarks of momentum p, the 
typical spatial separation | ai — y | is of order l/|p| (since x — y and p are Fourier 
conjugates in eq. (24)), and the closed loops that one would have to consider 
in the above argument have a typical area of order l/|p| 2 . On the other hand, 
numerical studies of the gauge fields produced in the McLerran-Venugopalan 
model indicate that the expectation value of Wilson loops decreases exponen¬ 
tially with the area of the loop when it becomes larger than Qj 2 [36, 37]. This 
suggests that this ambiguity should not affect much the quark spectrum for 

IpI ^ Qs■ 

8 This should be the case at least when x° —> +00 in a sensible model of the gauge fields 
produced in heavy ion collisions, thanks to the expansion and dilution of the system. 
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3 Statistical sampling method 

3.1 Sketch of a direct algorithm 

Eq. (22) contains the essence of our procedure for calculating the quark spec¬ 
trum: 


i. Draw randomly a pair of sources Jf 2 , and solve the classical Yang-Mills 
equation (3) with null retarded initial conditions. 

ii. For a given transverse momentum k±, wavenumber u, spin s and color a, 
initialize the spinor ip k± vsa as a free negative energy spinor. 

iii. Solve the reduced 2+1 dimensional Dirac equation (20) for the time evo¬ 
lution of this spinor over the color field found in the step i. 

iv. At some sufficiently large time, project the spinor 'ipk^wsa on a free positive 
energy spinor of transverse momentum p±, same wavenumber v 1 spin er 
and color b. This gives the reduced inner product [V’p+ I ,cr&|V’fc( Ll , S a] T ‘ 

v. Repeat the steps ii to iv in order to sum over all k±, v, s, o’s. 

vi. (optional) Repeat steps i to v in order to average over the configurations 
of the color sources Jf 2 . 

If we discretize the transverse plane into a N x x N y grid and the rapidity 
axis into a grid with N v points, the computational cost for solving the Dirac 
equation for one mode function scales as N x N y N T where N T is the number of 
timesteps. Note that this cost does not depend on N v since the 77 dependence 
can be factored out for a boost invariant background. This must be repeated 
for the N x N y N v mode functions. Therefore, the total computational cost scales 
as {N x N y ) 2 N v N T . 


3.2 Statistical method 


The algorithm described in the previous subsection is deterministic, but it suf¬ 
fers from an unfavorable scaling with the size of the transverse grid, since the 
computational cost scales as (N x N y ) 2 . Instead of doing in full the sum over all 
the modes, it is possible to do it by a Monte-Carlo sampling in which one uses 
random linear superpositions of the mode functions 


i’c 



d 2 k± dv 
(27t) 2 27T 


Ck±vsa 


$ 


k±vsa 5 


(28) 


where the coefficients c k± usa are random numbers with the following variance 
{ck ± vsaC* k ' ±v ' s 'a') = (27t) 3 5(k± - k' ± ) 6( V - v’) (5 SS / S aa / . (29) 
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The justification of this approach is to note that the projector on the subspace 
of negative energy spinors can be rewritten as follows 



d 2 k± dv 
(2ir) 2 2ir 


^kj_vsa) {^k^vsa 


J[Dc]P[c] |+)(+ 


(30) 


where V[c\ is a normalized probability distribution (in practice a Gaussian dis¬ 
tribution) with zero mean value and a variance given in eq. (29). In eq. (22), 
the integrals over k ±, v and the sums over s, a are then replaced by a statistical 
average over these random numbers, 


2tt<5(0) 


dN a 


1 


dy p d 2 p± 87 t( 27 t ) 3 


E 


=T,4. ' 

b 


dis 
27r 


[Dc]P[c] lim h |V> C ) 


(31) 

(L v is the total length of the 77 interval represented in the lattice implementa¬ 
tion.) Since the random sum in eq. (28) mixes 9 the various ids, the evolution 
of ip~ is governed by the 3+1 dimensional Dirac equation (16), and the com¬ 
putational cost of this approach scales as N x N y N v N T N con f where N con f is the 
number of samples used in the statistical average. Compared to the direct de¬ 
terministic method, a power of N x N y has been replaced by N con {, which is 
advantageous for large grids if N con f <C N x N y . 


3.3 Statistical errors 

The statistical method summarized by the eqs. (30) and (31) is exact only in 
the case of a perfect sampling of the Gaussian distribution V[c\. In practice, 
this sampling is performed by generating a large but finite number IV con f of 
configurations. In doing this, the left hand side of eq. (30) is replaced by 

^ E ^ 1 $}) (V>, | (32) 

J.J' 

where we have used discrete notations that correspond to the lattice implemen¬ 
tation, and we use the following shorthands : 

V = LxLyLr, (total lattice volume) 

J — {jxi jyi jrji d) • ( 33 ) 

(The integers j x ,y,ri label the momentum modes k_\_,v in the lattice implemen¬ 
tation, and s,a are the spin and color quantum numbers.) The coefficients 

9 One may replace eq. (28) by a random sum in which the modes v are not mixed. These 
restricted linear superpositions obey the reduced 2+1 dimensional Dirac equation (20), but 
one must repeat the resolution of the equation for each mode is, so that this modification has no 
merit in terms of computational cost. In fact, using eq. (28) and solving the 3+1 dimensional 
Dirac equation has the advantage that it is immediately generalizable to a non-boost invariant 
background color field if necessary. 
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that appear in the sum in eq. (32) result from N con f samplings of the Gaussian 
distribution : 

-i N conf 

c N „, (j, j 1 ) s e 4 I 4” • < 34 > 

iVconf t J J 
n —1 

(n) . -* 

In this equation, c4 is the n-th random sample for the mode J. 

With iV C onf samples, the observables we are interested in (e.g. the inclusive 
quark spectrum given by eq. (31)) are generically approximated by 

o Ncont ^ E > ( 35 ) 

/ j,j' 


where F denotes the quantum numbers of the final state, Ylf a partial sum 
over these quantum numbers (in the case of eq. (31), this partial sum is over the 
wave number v, the spin and color of the produced quark), and A f a normaliza¬ 
tion factor. C Nco M- J') is itself a random number, whose distribution can be 
determined in the large N con { limit by a method similar to the derivation of the 
central limit theorem. The mean value and variance of this approximation can 
be obtained from those of Cpf coa! (J, J 7 )- Let us first recall that 

(cf)=O t ( c fc^)=0. 

This leads easily to 

<C;w(J,j*)> = 

(C Nco M:J')CN cml (K,K')) = 


The formula for the variance is exact if the distribution of cj is Gaussian. Like 
in the central limit theorem, the variance of Cjv conf decreases as l/lV con f. 

From the first of eqs. (37), we obtain the mean value of Ojv conf 

<»«-> =*£ WWW) • < 38 > 

/ J 


(cf*c { p) = V6 nn ,6 f f , . (36) 


V i 

(C Nco M,f))(C Ncon{ (K,K')) 


(37) 


which is indeed the exact value of the observable. Using the second of eqs. (37), 
we get 


( 0 N C onf) 


(<^ conf ) 2 

+ n^v^ E E W£) ■ 

/,/' J K 

(39) 


We see that the standard deviation of On coo! decreases as l/y/N con {, with a 
coefficient that has a non-trivial covariance. It can itself be estimated by the 
statistical method as follows : 
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i. Define two random linear superpositions of the negative energy mode func¬ 
tions : 

^2 ( 4 °) 

J 

with uncorrelated random weights Cj and Cj\ 

ii. Evolve these two spinors in time by solving the Dirac equation, 

iii. Compute the following quantity : 

N 2 (^? + 1) (^2 1 i ’ 0 /) f 2 HvjfVV) = 

= (V’ 2 _ l^ + )(^ + |^r)| 2 , (4i) 

/ 

iv. Repeat the steps i-iii in order to average over the random numbers Cj . 
Since this is just an error estimate, a small number of samples is sufficient. 
In practice, one may divide the N con { samples already calculated in two 
subsets, and use these subsets to evaluate the error. 

The standard deviation of C ) 7 v conf is the square root of the result of this com¬ 
putation, divided by \JN con {. Note that the summand in eq. (41) is a complex 
number, which can lead to phase cancellations when summing over the final 
quantum numbers /. These cancellations are more effective for more inclusive 
observables, thanks to a more extended sum on /. 

3.4 Relation to “low-cost fermions” 

In real-time lattice simulations of fermions, the so-called low-cost fermion me¬ 
thod [38] has been used in several works, e.g. [39, 40, 41, 42, 43, 44], Let us 
briefly compare this method with our approach. In our statistical method, we 
compute the following quantity using the stochastic field (28): 

(C f (x)OC(!/))c = ^2^}\x)Oipy{y) , (42) 

J 

where J comprises all quantum numbers including momentum, and O is a 
matrix that depends on the observable we wish to compute. In the case of the 
spectrum (31), O = (x)ip 0 ^ \y). In the low-cost fermion method, instead of 

using one stochastic field (28), one employs two kinds of stochastic fields called 
“male” and “female” fields: 

Vv = 53 [ c j^/ + dj^s] - Vv = 12 l c J^’S - d J^~j\ ’ ( 43 ) 
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where c j and dj are independent random numbers that have the same variance 
as eq. (36). Combining these two Helds, one can compute 


-(vifaWp ( y )> 




instead of eq. (42). By using the completeness relation 


^/^>x)^ f (Cy) +V' J -(Cx)^ / t (t, y) 


= 1 


(44) 


(45) 


(1 is the unit matrix in the spin, color and position of the spinors at the time <), 
we can relate the quantities evaluated in our method (42) and in the low-cost 
fermion method (44) by 

(C t ( a: ) 0 C(y)} = -(4( I ) 0 iW) + 2i ( x “y) • ( 46 ) 

(for 2 spin states and 2 colors.) Therefore, the two methods provide the same 
result 10 up to a trivial additive term, that can be interpreted as a constant vac¬ 
uum contribution. However, our method has two advantages over the low-cost 
fermion method. Firstly, it is numerically less costly than the low-cost fermion 
method, simply because it uses only one kind of stochastic field. Secondly, the 
statistical errors are smaller for the evaluation of the spectrum. In our method 
(42), the spectrum is directly obtained from the statistical ensemble, without 
any subtraction. On the other hand, in the low-cost fermion method (44), one 
gets directly access to | — / (/ being the fermion occupation number), and the 
vacuum 1/2 must be subtracted later. Because this vacuum 1/2 also contains 
statistical errors due to the Monte-Carlo sampling, the low-cost fermion method 
suffers from comparatively larger statistical errors, especially when the value of 
the occupation number is small compared to 1/2. 


4 Fermionic mode functions on the light-cone 

4.1 Background gauge field 

In order to use the classical-statistical method in the calculation of the quark 
spectrum, we should solve the Dirac equation with a background color field, 
starting with a free spinor at = — oo. However, it is not straightforward to 
do this numerically, due to the singular nature of the gauge field on the light- 
cones x ± = Q. The field strength on these lines is proportional to a d)^), 
which cannot be handled easily in a numerical program. Similarly to the case 
of gluon production, one should first solve the Dirac equation analytically up 
to a surface r = To , and perform the numerical resolution only in the 

forward light-cone for t > t$. 

10 This is not the case if the initial state is not charge neutral. 
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In this section, we generalize to the case of spinors the derivation that was 
performed in ref. [45]. Since the Dirac equation is linear, its solution can be 
written as the sum of a left-moving and a right-moving partial waves, as illus¬ 
trated in the left diagram in the figure 1. These two partial waves are totally 
independent. We take advantage of this fact to choose the gauge for the back¬ 
ground field differently for each of them, in order to simplify the resolution of 
the Dirac equation. Of course, before adding up the two partial waves in order 
to construct the full solution, we must perform a gauge rotation that brings 
them to a common gauge. 



Figure 1: Left : decomposition of the solution of the Dirac equation into left- 
moving and right-moving partial waves. Right : structure of the gauge field 
produced by the two nuclei in the A~ = 0 gauge, that we use for solving 
the Dirac equation for the right-moving partial wave. denote light-cone 
coordinates, x^ 1 = (t ± z)/V2. 

For the right-moving partial wave, it is convenient to work in the light-cone 
gauge A~ = 0, in the same spirit as what was done in ref. [45] for the gluons. 
After having determined the right-moving spinor at the proper time r 0 inside 
the forward light-cone, we will rotate it to the Fock-Schwinger gauge. The 
calculation for the left-moving spinor (to be performed in the light-cone gauge 
A + = 0) will not be performed here, since the result can be guessed from the 
other partial wave by symmetry. The structure of the background field in the 
A~ = 0 is shown in the right diagram of the figure 1 . It is made of the following 
two elements : 

1. the nucleus moving in the +z direction produces a field A ^, proportional 
to a S(x~) and independent of x + . 

2. in the region x + > 0,x~ < 0, the nucleus moving in the —2 direction 
produces a field A\, which has the form of a transverse pure gauge, 

A. = - g U\d l U 2 . (47) 
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3. in the strip corresponding to the shock-wave of the nucleus moving in the 
+z direction, this field A\ receives a color precession induced by the first 
nucleus. This color rotated field reads 

< = l -U lab (x- lXl _) {U\d l U 2 ) b , (48) 


where 


U\{x ,jc_l) = T exp ig / dz Af(z ,x±) . 
Jo 


Note that eq. (48) is equivalent to 


4 = = UiA^ul . 


(49) 


(50) 


Therefore, by starting the evolution at x° = — oo, the right-moving partial wave 
encounters first the field A 2 and then the held A+. In order to express the quark 
spectrum, we need the negative energy spinors, r/v (fc is the 3-momentum of 
the incoming quark, s its spin and a its color). Before they encounter the 
background color held, they read simply 

= (51) 

Since the background color held has only a hnite jump on the half-line x + = 
0, x~ < 0, the spinors are continuous across this line, and the above formulas 
remain valid up to x + = 0 + (just above this line). 


4.2 Evolution in the region x < 0,a: + > 0 

The next step is to solve the Dirac equation in the region x + > 0, x~ < 0. Since 
the background held is a pure gauge in this region, the covariant derivative can 
be written as 

= U\ d» U 2 . (52) 

Therefore, the new spinor dehned by tpkso, = ^2 V 1 ksa °b e y s the free Dirac equa¬ 
tion : 

{i0 - m )i>ksa = 0 • ( 53 ) 

The solution of this equation can be expressed in terms of the initial value of 
V’fcso on the surface x + = 0 + by the following Green’s formula, 


^ksa( X ) 



dy d 2 y ± S° R {x,y)j+-ij ksa (y) , 


where 5° (x, y) is the bare retarded propagator for a quark, 



d 4 p g-ip-Oc-I/) 
(27t) 4 ]/ — m + fpo7°e 


(54) 


(55) 
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Note that, even if we have not restricted the integration range for the variable 
y~ in eq. (54), the fact that the support of the retarded propagator 5° is limited 
to £° > y°, (x — y) 2 > 0 imposes that y~ < x~. Therefore, for a point x located 
in the region x + > 0 , x~ <0 (shaded in green in the figure 1 ), only points with 
y~ < 0 on the initial surface can contribute. This is the reason we can solve 
independently the equation for the left- and right-moving partial waves. 

By introducing the Fourier representation (55) of the retarded propagator in 
the Green’s formula (54), one can perform most of the integrations analytically 
except the final integration over transverse momentum. This leads to 


= uUxx)J^e ip ^e^ +x - + S^U 2 (p^+k ± ) 

(. + p*7 l + r?A 


1-7' 


2 k+ 


V + v s (k) 


(56) 


where M p = \J + m 2 is the transverse mass and where we have introduced 
the projector V + = 7 _ 7 + / 2 . We have also introduced the Fourier transform of 
the Wilson line U 2 


U 2 (k ± ) = j d 2 y± e" 4 *™ U 2 (y±) . 


(57) 


One can check that eq. (56) falls back to the free spinor v a (k) exp(ifc • x) if the 
background field is turned off (U 2 = 1). 


4.3 Evolution across the line x =0 ( x + > 0) 

The next step is to propagate the spinor (56) across the field Af of the nucleus 
moving in the +z direction. The region of space-time supporting this field is 
infinitesimal (since A + ~ <5(rr — )), but the infinite strength of the gauge field in 
this region nevertheless produces a finite change of the spinors. In this region, 
there is also an 0(1) transverse component of the gauge field, given in eq. (48). 
The Dirac equation, 

[?'(<9 + — igAf)"/~ + id~ 7 + — iD l 7 * — m] if>~ = 0 , (58) 

can be separated into a part that depends on the background field A+ and 
a constraint independent of Af by multiplying 11 it by the projectors V + or 
V~ =7 + 7“/2: 


id ^ + ^ksa 


i{d + - igA+)V ip ks , 


m-vfD 1 _ 

- 2 - 'y V ksa 

7 + (*7*D I +m) 
- 5 - Vksa ■ 


(59) 


The first equation, independent of the background field, is a constraint that 
relates the two projections of the spinors at every x~ (note that this equation 

11 One may use the following identities : = V~ 7 — = 0, V~ 7" 1 = 7' 1 ^ and V +7 = 

7 ~V~. 
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does not contain the d + derivative). The second equation determines the dy¬ 
namical evolution (in the variable x~ ) of the V~ projection under the influence 
of the background field Af. Inserting the first equation into the second gives a 
second order equation that drives the evolution of the V~ projection, 

(2d~ (<9+ - igA+ )-D 2 ± + m 2 ) V~ip ksa = 0 • (60) 

In this equation, the d + derivative and the field A ^ are large (inversely pro¬ 
portional to the thickness of the shock-wave that supports the A+ field), while 
all the other terms do not have this large factor. Physically, keeping only the 
term in d + — igA + leads to the eikonal approximation, where the fermion would 
propagate on a straight line along the x~ axis, while the terms — D\ + m 2 lead 
to some transverse diffusion with respect to this axis. In the limit where the 
thickness of the shock-wave goes to zero, we can neglect it : 

i(d + - igAt)V~ip ksa = 0 . (61) 

(Note that this approximation is only valid to cross the shock-wave, and should 
not be used to evolve at a finite distance from the shock-wave.) The solution of 
this equation is very simple, 

V ~^ksa( x ) = Ui{x-,x ± )V~^ sa {0,x + ,x ± ) , (62) 


where the ^“-dependent Wilson line U\ was defined in eq. (49). The spinor 
at x~ = 0 that appears in the right hand side is given by eq. (56). Then, the 
constraint equation can be solved to give 


VA’ksaix) 



m — i'fD 1 
2 


7 U x {x ,x ± )V ip ksa (0,z + ,x ± ) 


(63) 


Note that the constraint defines the V + projection of the spinor only up to an 
arbitrary function of x~ and This “integration constant” can be determined 
by requesting that we recover the V + projection of a free spinor when all the 
Wilson lines are set to the identity. 


4.4 Transformation to Fock-Schwinger gauge 

In order to gauge transform these spinors to the Fock-Schwinger gauge, it is 
sufficient 12 to multiply them by U\, 

VW FS (*) = U \ f(*-L)^,a(*) • ( 64 ) 

When this transformation is applied to eqs. (62) and (63), the Wilson line U\ 
appears in two types of combinations : 

U\U 1 = l and U\D i U\ . (65) 

12 In the case of the background gluon field, an additional transformation was necessary, see 
the eq. (18) of ref. [45]. But since the color rotation that characterizes this transformation 
is equal to the identity for proper times r <C Q7 l (see the eq. (19) in ref. [45]), it has no effect 
on the spinors. 
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The second of these structures can be simplified if we recall that D l is the 
covariant derivative built with the field of eq. (50). Therefore, 

U\D i U l = U\ (<9* — igUiA l 2 Ul)Ui 

= d i - ig(A\ + 4) , (66) 

s» _✓ 


where A\ is defined in the same way as A\ (see the eq. (47)), 


A[ = l - U\&U i . 


(67) 


Note that the field A\ + A\ that appears in this equation is nothing but the 
transverse component of the gauge potential in the Fock-Schwinger gauge at 
t = 0 + , hence the notation D* for the resulting covariant derivative. Therefore, 
the two projections of the right-moving negative energy spinors in the Fock- 
Schwinger gauge read (at x~ = 0 + , just above the shock-wave) 


'P Vw. P s(*) = U k x ±) 


d 2 p± 


(2tt)2 

v+ ^ksa FS ( x ) = (^-Dps - m) U\ (x ± ) 


e' p e'asS 

p l 7 * — m i 

_ r\i' 

d 2 p± 


x U 2 (p± + k ± ) - ^ fc +"‘ 7 +u s (fc) 


(2tt) 2 


1 p ^ x -l 


x ^(pj. + k l)7 


_ p 7 — m , 


2 M 2 


7 + u s (fc) (68) 


The eqs. ( 68 ) for the right-moving spinors must be completed by a set of similar 
equations for the left-moving spinors. These can be obtained from the above 
formulas by the following substitutions 


U 2 ->• Ux 

x + —> x~ , A; + —>■ fc _ 

v + ^ v~ , v~ ^ v + 

7 + —> 7 ^ , 7~ —> 7 + . (69) 


At this point, we have the components of the negative energy mode functions on 
the light-cone r = 0 + (i.e. just after the collision), which provides all the nec¬ 
essary initial data for studying their evolution after the collision. As explained 
in the appendix C, these formulas can also be used as initial conditions at some 
initial time r 0 > 0 , provided that r 0 a± where a± is the transverse lattice 
spacing used in the numerical resolution. 


4.5 Mode functions in the v basis 

So far, we have derived the mode functions i^ksa in terms of the Cartesian 3- 
momentum k. However, as explained in the section 2.3,the boost invariance of 
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a high energy collision is more manifest if one uses the modified spinors defined 
in eq. (13) and if one further goes to a basis where the spinors are labeled by the 
quantum number v (Fourier conjugate to 77 ) instead of y. It is easy to obtain 
these new mode functions by the following transformation : 

/ + OO 

dy e lvy ip k±vsa ■ (70) 

-OO 

For the right-moving partial waves, the integral that enters in the transformation 
of V e il>- is of the form 


(■ + 00 




dy e wy e £ 5 e %ae 
while for the left moving partial waves, one needs 


r+00 


r = 


dy e wy e e i e iae 


(71) 


(72) 


These integrals can be expressed in terms of the F function, 

P R = {-iaT + i r(-«/-§) , I[ = T(iv + |) . (73) 

Let us now recapitulate our results for the fermionic mode functions on the 
light-cone after this transformation, after summing the right-moving and left- 
moving partial waves, and including both the V + and V~ projections 

{ / jw -2 \ iis + 

r (-^+g) c/ 2 t ( ;K J-)C72(p_L + fej.)7 + 

+ e_ ^ r {2^) r ( iu +l) u l( x ±) u i(p±+ k ±)'y~np l 'y l + m )v s (k Xi y=o). 

(74) 


In this formula, we have kept only the terms that are non-vanishing in the limit 
t —» 0 + . Therefore, its use should be restricted to very early times. 


5 Summary and outlook 

In this paper, we have reconsidered the problem of quark production in heavy 
ion collisions in the color glass condensate framework. Our approach is closely 
following that of refs. [31, 32], which expresses the leading order inclusive quark 
spectrum in terms of a set of mode functions of the Dirac equation, but over¬ 
comes a number of limitations of this earlier work. 
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Firstly, we have defined a basis of fermionic mode functions that are more 
appropriate for the boost invariant expanding geometry of a high energy col¬ 
lision. In particular, since these mode functions are indexed by the Fourier 
conjugate v to rapidity, they are especially suitable for a lattice implementation 
in which one discretizes the rapidity axis. We have calculated analytically the 
value of these mode functions just after the collision, at a proper time Q s t <C 1 
and in the Fock-Schwinger gauge, in terms of the Wilson lines that represent 
the classical color background field of the two colliding nuclei. Thanks to these 
analytical initial values, one will not have to deal with crossing the light-cones 
in the numerical resolution of the Dirac equation. 

Secondly, we have exposed a statistical method for sampling the set of modes 
over which one must sum in the calculation of the quark spectrum. This method 
ensures that no mode is left out, while considerably reducing the computing time 
compared to a complete sum over all the modes. This approach also provides a 
more robust way of estimating the error one makes in the sum over the modes 
functions. 

In a forthcoming paper, we will apply the formalism that we have setup in the 
present paper to a study of quark production in two situations. Firstly, we will 
present a test of the method in the case of a background field for which one can 
solve analytically the Dirac equation for the mode functions (a constant SU (2) 
chromo-electrical field). Then, we will present results on quark production in 
heavy ion collisions, in the case where the background color field is given by 
the MV model. In order to mitigate the problems caused by the lattice fermion 
doubler modes, we will also explain how our framework must be modified in 
order to include a Wilson term in the fermionic action. 
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A Quark spectrum in the Schwinger-Keldysh for¬ 
malism 

In this appendix, we recall the derivation of the formulas (4) for the inclusive 
quark spectrum, starting from the standard LSZ reduction formulas. For sim¬ 
plicity, we consider only one flavor of quarks. Since in the CGC framework, the 
external sources are only coupled to gluons and quarks can only produced in 
quark-antiquark pairs (the net flavor number is always zero). 
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A.l Single quark pair amplitude 

Firstly, let us consider the amplitude for producing a single quark-antiquark 
pair, 

= (Q(p),Q(q) out|O in ) 

= (dout | ^out (p)^out ( 9 ) |0in) • (75) 

In order to keep the notations concise, we are not writing explicitly the color and 
spin indices of the quark and antiquark. The operator bl ut (p) (resp. d' out (q)) 
creates a quark of momentum p (resp. an antiquark of momentum q). 

Using the decomposition of the free field operator ipout(x) as a superposition 


of free modes, we have 



&out (p) = 

J d 3 x u(p) 7 0 tfrout(U x) e lp ' x 


dout(q) = 

j d 3 x ^ out (t,x) 7 0 v(q) e lq x . 

(76) 


(In these formulas, p° = \J pi 2 + m 2 and q° = \J q 2 + m 2 .) Note that the time 
t at which these formulas are evaluated do not change the result. Using these 
formulas, standard manipulations lead to the LSZ reduction formulas for the 
single quark pair production amplitude, 

Mi(p,q) = 

x(0 O ut|TV’(ar)^(j/)|0 in ) (i 0 y +m) v{q)e iq ' v . (77) 

Note that the 2-point correlation function that appears in this formula is a 
Feynman (i.e. time-ordered) propagator. Its evaluation to all orders in the 
background gluon field cannot be performed in practice. Indeed, even though it 
obeys the Dirac equation, its determination is made extremely complicated by 
the fact that it must satisfy mixed boundary conditions, both at x° = —00 and 
at x° = +00. 

A.2 Inclusive quark spectrum 

There is no practical way to calculate the single quark amplitude considered in 
the previous subsection, in the presence of a strong background gluon field, as 
is the case in the high energy heavy ion collisions. This is not a big limitation 
however, because this quantity is also not very phenomenologically useful in 
such a context. Indeed, for light quarks (quark flavors for which m < Q s ), quark 
production is not a rare process and more than one quark pair are produced in 
a collision. Therefore, the probability Pi of producing exactly one quark pair 
(given by the square of M. 1 ) is not very interesting. 

Much more useful would be the complete probability distribution, Pi, Pi, P 3 , 
etc. Unfortunately, calculating them is as complicated as calculating P\. But it 
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turns out that the moments of the probability distribution are much easier to 
compute. Besides the trivial one (X)(^Lo -Fn = 1), the simplest of these moments 
is the first moment, Xm n Pt; that gives the mean number of produced pairs. 
A little more information can be gathered by considering the same quantity in 
differential form, which is precisely the quark spectrum. In terms of transition 
amplitudes, this quantity reads 


dN q 

d 3 p 


1 +oo 1 r n 

5 <n + J 

X | (Oout |^out(p)d out (q)b out (Pi)d out (<Zi) ' ^out (Pn)d out (<?n)|0 in )| 


(78) 


where we have used the shorthand 


d$ q = 


d 3 q 

(2tt) 3 20J q 


(79) 


for the invariant phase-space of final state quarks and antiquarks. In this for¬ 
mula, the differential probability for producing n + 1 quark-antiquark pairs is 
weighted by the number of quarks (n+ 1), and then integrated over the phase- 
space of all the antiquarks and of n of the quarks. This quantity is normalized in 
such a way that it gives the mean number of produced quarks after integration 
over d 3 p, 

J d3p = £ nPn - ( Nq) ■ (so) 

Most of eq. (78) is in fact the projector on the subspace of states with net flavor 
number — 1, 


+oo i /» n 

n—0 v ' J i= 1 

x | [pi p n ] Q [qqi ■ ■ ■ q,<_ 0 ou t)([pi • • • Pn] Q [qqi ■ ■ ■ q n ]q out | , (81) 

that has a trivial action on the state &out(p)|0; n ) (since this state has also flavor 
number —1), and eq. (78) can then be reduced to the much more compact form 


d 3 p ~ (2ir) 3 2oj p <M & °ut(P)Wp)|0i„> • (82) 

The interpretation of this formula is quite transparent, since it amounts to 
evaluating the expectation value of the final quark number operator, for a system 
prepared in the pure state |0i„). 

Using eqs. (76), this can be rewritten in terms of the quark field operator as 
follows, 

S = Jd'^ve^uip) (i ?,-m) 

x(O i „|V>(a;)'0(j/)|Oin) {i $ y -m) e~ ip v u(p) . (83) 
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The main differences between this expression and the similar expression for the 
amplitude A4i are the following : 

1. The vacuum state is the in-vacuum state on both sides. 

2. The two spinors are not time ordered. 

These differences in fact lead to considerable simplifications in the evaluation of 
this quantity in the presence of a strong background gluon Held. 

The first step is to note that the 2-point correlator that appears in the 

integrand of eq. (83) is the component S _(_ (ar, y ) of the fermion 2-point function 

in the Schwinger-Keldysh formalism. In the presence of a background field A p , 
its tree level expression (to all orders in the background field) can be obtained 
by noticing that it obeys the following equations 

[i t> x —m)S — \-(x,y) = 0 , S^+(x,y) [i t> y -m) = 0 

lim S„ + (x,y)=S v —(x,y) . (84) 

x u ,y' J —> — oo 

Next, one should recall the expression of the vacuum propagator S^ uunl (x,y), 

Sl—(x,y) = E / (2^)32^ etHX ~ V) Mk)v s (k) ■ (85) 

It is then easy to construct a senri-explicit expression for the dressed propagator 
S _ \-{x,y) in terms of a basis of solutions of the Dirac equation: 

/ r fiu _ 

j2n^J’ ksa{x) ^ a {y) 

(i0 x -m)ilj ksa (x) = O , lim ip ksa {x) = v s (k) e lk ' x . (86) 

x°—> — oo 

When this expression is inserted into eq. (83), one must evaluate the following 
expression 

j d 4 x e zp ' x Uo(jp) {i 0 X -m) ^ ksa {x) . (87) 

It is useful to note that 

e ip x u<j(p) (i $ x +to) = 0 . (88) 

Adding this identity to the integrand of eq. (87), we obtain 

J d 4 x e ip x u a {p) (i $ x to) ip ksa (x) = i J d 4 x d^ [e ip ' x Uaiph 11 i’ksaix)] , 

, (89) 

which can be rewritten as a 3-dinrensional integral since it is the integral of a 
total derivative. The boundary at spatial infinity can be dropped if there are 
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no background fields there. The boundary at —> —oo does not contribute 
because it leads to the vanishing overlap v) (p)v(p). The only remaining contri¬ 
bution comes from x° —> oo, 


[ d 4 x e lpx u a (p) (i 0 X - m) ipksa(x) = i lim 
J x °—>-+oo 

This formula leads immediately to eq. (4). 


d 3 x e lpx u\{p) tp ksa {x) . 

(90) 


B Quark spectrum from Feynman amplitudes 

In the previous appendix, we presented a derivation of the quark spectrum based 
on fairly standard many-body manipulations: we first related this spectrum to 
the expectation value of the quark number operator, and then we evaluated 
the latter using the Schwinger-Keldysli formalism. However, a more elementary 
derivation is possible, where the many-body aspects of the problem are treated 
“by hand”. In the present appendix, we present such an alternate derivation, 
starting from the Feynman amplitudes for producing 1,2,3,... quark-antiquark 
pairs, and combining them in the appropriate way to obtain the single quark 
spectrum. This method is a bit more involved since it requires to account for all 
the final state particle permutations, but it has the advantage of making more 
tangible the combinatorics that happens under the hood in the derivation of the 
appendix A. 

B.l Pair production amplitudes 

The starting point is the amplitude A4i(p, q) for producing one quark-antiquark 
pair, already introduced in eq. (75). This amplitude is made of a time ordered 
2-point function connecting the quark of momentum p and the antiquark of 
momentum q , times a disconnected sum of vacuum-vacuum graphs. The latter 
is crucial in the presence of a background field, since the sum of the vacuum- 
vacuum graphs is not a pure phase (unlike when the background is the vacuum). 
In practice, we do not need to calculate this factor, since it is the same in all 
amplitudes and can therefore be determined at the very end by the request that 
the sum of all probabilities to produce 0,1,2,3,... quarks is equal to one. For 
now, we will simply write 

Mi{p,q)= (0 ou t|0in) x Ml{p,q) , (91) 

sum of vacuum graphs 


where is the connected part of the pair production amplitude (only this 
factor carries a dependence on the momenta of the produced quark and anti¬ 
quark) . 

In the rest of this appendix, we limit the discussion to the lowest order 
for the factor At this order, it is simply made of a Feynman propagator 
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connecting the produced quark and antiquark, dressed by the background field: 


M c 1 (p,q) = 



(92) 


For the sake of simplicity, we will represent this dressed propagator as follows: 



A more explicit expression of the connected part of the single pair production 
amplitude is 


M\(p,q) = u(p)T F (p, -q)v(q) , (94) 

where T F is the dressed Feynman propagator, amputated of its final free prop¬ 
agators. In terms of the dressed (G F ) and free (G°) Feynman propagators, it 
is given by 

G f =G° f +G° f T f G° f . (95) 

The convention for the momenta in eq. (94) is that the 4-momentum — q enters 
on one side of the propagator and the 4-momentum p exits on the other side. 

At the same level of approximation, the amplitude for producing n pairs is 
obtained as follows: 

M n (pi • ■ Pn, qi - q n ) = (0 ou t|0in) ^ e(cr) M\{p\, q ai ) • • • Mi(p n , qa n ) , 

cr£©n 

( 98 ) 

where & n is the symmetry group of the set [l,n]. The sum over all the per¬ 
mutations a £ & n is necessary in order to account of all the possible ways to 
connect the quarks with antiquarks. e(<r) is the signature of the permutation a, 
resulting from the signs collected when permuting fermion fields. 

B.2 Final state combinatorics 

It is possible to encapsulate all the information about the distribution of the 
produced quarks in the following generating functional, 

OO 

F[z(p)\ = ^2 J Z ^' Pl \u Z ^ Pn ' > \ M n{Pl---Pn,qi---q n )\ 2 , (97) 

Tl—0 ... p n 

<n -qn 
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where we have used the following shorthand for 1-particle phase-space integrals: 


+ 



p 


d 3 p 

(2tt) 3 2ui p 


d A p 

(2 


2i t9(p°)S(p 2 


TO 2 ) . 


(98) 


The + superscript on the integration symbol indicates that we keep only the 
positive on-shell energy. Likewise, a — superscript will indicate that the negative 
on-shell energy is retained: 


/ =J 2 tt d(-p°)S(p 2 - to 2 ) . 


(99) 


P 


From its definition, it is easy to see that the single quark spectrum is obtained 
as 


dN q 5T[z\ 
d 3 p 5z(p ) Z = 1 


( 100 ) 


Note also that unitary (the sum that all probabilities should be one) implies 
that T[z = 1] = 1. Inserting eq. (96) into eq. (97), we obtain 


+ 

?[ z \ = |(0out|0in)| 2 y] ^2 J Z (Pl) ■ ■ ■ Z (.Pn) 

n a,a'e.& n p 1 -p n 

<11 ■ -Qn 

xMl*{p 1 ,q ai )Ml{p 1 ,qy i ) ■ ■ ■ Mf(p n , q„ n ) M^(p n , q a > n ) (101) 

When all the spin and Dirac indices are summed over, the product in the second 
line forms closed quark loops, from which the spinors can be eliminated by using 
u(p)u(p) =jJ+m and v(q)v(q) = { — m. It is convenient to change q —> —q for 
all the antiquarks, so that we have: 




|<0ou t |0 in )| 2 £ 


(-l) r 


pee„ 


YlHz] 


QiQpi 


Qi Qr 


i =1 


( 102 ) 


where we have defined 13 


+ 

L[z]qq' = J z ( p ) T p(?,p)^+m)T F (p,g , )(/ + m) = 

p 



In the diagrammatic representation used for this quantity, the dotted line rep¬ 
resents the final state. Right of this line is an amplitude and left of this line is 

13 We denote T* ( q , p) = 7 0 Tj (p, q ) 7 0 . 
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a complex conjugated amplitude. In order to go from eq. (101) to eq. (102), we 
have used the symmetry of the n-quark and n-antiquark phase-spaces and the 
permutation p is defined as p = <j~ l cr'. 

Every permutation p £ G n has a unique decomposition in a product of 
cycles 14 . In eq. (102), each of these cycles will produce a closed quark loop, 
that depend on 2 through the quantity L[z] (linear in z) defined in eq. (103). 
The degree in z of such a quark loop is the order of the cycle (the number of 
iterations before the cycle returns to the starting point). For a cycle of order r, 
we will denote the value of the corresponding quark loop tr ((L[z]) r ), where the 
trace symbol compactly encapsulates the integrals over all the momenta along 
the loop, as well as the contractions over Dirac and color indices. The important 
point is that in eq. (102) the product of L’s depends only on the orders of the 
cycles into which the permutation p can be decomposed: if p = c\C 2 ■ ■ ■ Ci where 
the Cj ’s are cycles of orders rq, • • • ri , then we have 




j=i 


(104) 


In order to perform the sum over the permutations p, it is sufficient to know the 
number of p’s that admit a decomposition into aq cycles of order 1, a -2 cycles of 
order 2, ..., a n cycles of order n (with the constraint ai + 2 a 2 H-+ na n = n), 


n\ 1 

ad • • • a n \ l°i • • • n an ’ 


(105) 


and its signature 

n 

e(p) = (-i)"II(- 1 ) 0i • ( 106 ) 

1=1 

14 We recall the reader that a cycle is a circular permutation 1 —» —> (rra) i 


28 



Combining these results, we can rewrite 15 


E^E'M / fl^w.= 


n>0 


pee n 


Ql“-Qr 


i=l 


=e e 


— ao! 

n>0 a 1 +2a 2 + --- j = 1 J 


= E E 


i> 0 ai+ 2 a 2 H—=n j=l 

-e e n^(- 

p>0 0l+02+-'=P j = l J \ 


=E^ |-‘ r E 


(W 


p> 0 


tr (my 


= exp (tr ln(l — L[z])) . (107) 


j =i 


The second and third lines are equivalent because the constraint 5Z ? > i 3 a j = n 
prevents a^’s with j > n from being nonzero. Going from the third to the fourth 
line merely corresponds to a different way of slicing the sum over all (ij 's. 

The only missing ingredient is the prefactor |(0 O ut|0in)|C that can be triv¬ 
ially determined in order to satisfy unitarity. The final expression for the gen¬ 
erating functional is therefore 


exp (tr ln(l — L[z])) 
exp (tr ln(l — L[l])) 


(108) 


Taking a functional derivative leads to the following expression for the quark 
spectrum: 


dN q 

d 3 p 



- CD 


At this stage, we have a representation of the quark spectrum in terms of the 
object L[z\, which is itself quadratic in the Feynman propagator of a quark 

15 The final expression is reminiscent of the determinant of a Dirac operator, and could 
probably be derived more straightforwardly by path integral methods. 
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in a background field. However, this expression contains terms of arbitrary 
order in this propagator, because of the prefactor (1 — L[l]) _1 . Note that in a 
weak background field, the quantity L[ 1] is small and the first term, quadratic 
in the Feynman propagator, dominates the sum. This is the limit where the 
pair production probability is small and where at most one pair is produced in 
a collision. Therefore, the quark spectrum is almost equal to the differential 
probability of producing a single quark, given by the first diagram in the above 
series. The second, third, etc... diagrams encode Fermi-Dirac correlations that 
are only important when more than one quark are likely to be produced (which 
is the case in heavy ion collisions for quarks whose mass is comparable to the 
gluon saturation momentum or smaller). 

B.3 Expression in terms of retarded amplitudes 

It turns out that a much simpler expression, quadratic in the propagator, can be 
obtained if we rewrite eq. (109) in terms of the retarded propagator of the quark. 
One can define a retarded analogue T R of T F , built with free retarded propa¬ 
gators instead of free Feynman propagators. The free Feynman and retarded 
propagators are related by 

G°(p) = G°(p) + 2Tr(tf+m)9(-p°)5(p 2 -to 2 ) . (110) 

'-V-' 

p-(p) 

If we denote by V one insertion of the background field, the following two 
equations define T F and T R recursively 

T f = V + VG° f T f = V + T f G° f V 
T r = V + VG° R T R = V + T R G° R V , (111) 

from which one first obtains: 

t f = (i -vGiyW 

= (1 - vg° f )-\ 1 - vg° r )t r = (1 - vg° r Vp_)-\ 1 - vg° r )t r 

= (1- (1 -VG° R )- 1 Vp_y 1 T R = (l-T R p_)- 1 T R . (112) 

In terms of these compact notations, L[ 1] can be written as 

m = T* P+ T fP - (113) 

where p+(p) = 27 + m)9(+p°)S(p 2 — m 2 ), and manipulations similar to above 
lead to an expression that depends only on the retarded T R : 

L[l} = (l-T* R pJ)-\T* R p + T R pJ)(l-T R pJ)- 1 . (114) 

Note that since L[z\ is a linear functional of z(p), the derivative 5L[z]/5z(p) is 
similar to L[ 1], but with the intermediate momentum p fixed instead of being 
integrated over. We will denote this as follows: 

= (1 - T* rP -)-\t* r (rf+ to) T rP _) (1 - T.p.y 1 . (115) 
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In the same fashion, one also obtains the following identities: 


(l — T rP~) (l — TrP-) 
1 - L\ 1] 


1 + (T*p + T R p_) , 


Therefore, 


dNq 

d 3 p 


-tr 



6L[zy\ 

Sz(p)J 


-tv{T* R {j/+m)T R p_) . 


(116) 


(117) 


The outcome of these algebraic manipulations can be pictorially summarized as 





(118) 


where the blue lines in the right hand side are retarded propagators. The con¬ 
nection with the previous appendix and the formula (4) resides in the following 
identities: 


-tr(r*G/+m)T>_) 


Ua(p)T R (p, -q)v s {q) 


f K ( p ) t «( p >- 9 K (<?)| 2 , 

f7,S 

lim [ d 3 x x) x) . (119) 


C Evolution from r = 0 + to To > 0 

The eqs. (68) (and their counterparts for the left-moving partial waves) give 
the value of the fermionic mode functions immediately after the collision, one 
the semi-axis x~ = 0 + ,cc + > 0 for the right-moving partial waves and on the 
semi-axis x + = 0 + ,aT > 0 for the left-moving partial waves. Therefore these 
formulas provide for each mode function the complete initial data on the light- 
cone r = 0 + . 

However, the transformation to the quantum number v (Fourier conjugate 
of the spatial rapidity rf) introduces a non-analyticity in the time dependence, in 
the form of factors t ±w . For this reason, the numerical resolution of the Dirac 
equation should be started at a strictly positive proper time To > 0. Therefore, 
one should evolve the mode functions from the time r = 0 + to the time tq (see 
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Figure 2: Evolution between the forward light-cone (r = 0 + ) and the time tq 
at which the numerical solution of the Dirac equation starts. 


the figure 2), by using the Green’s formula for solutions of the Dirac equation, 


fp{r 0 ,ri,x ± ) 


1 d, * + 1 
'-y— —0+ y+=0 + 

!/+>0 y ~> 0 



d 2 y _l S r (tq, 77 , x±;y) iji i/j(y) . 


( 120 ) 

This Green’s formula in principle contains the Dirac propagator dressed by 
the background color field, which is not known analytically inside the forward 
light-cone. However, if one is interested only in very small propagation times 
Q s to -C 1 one can use the bare Dirac propagator instead, the effects of the 
background field on the evolution of the spinors start becoming important only 
at Q s t 0 > 1. But even with a bare propagator, the evaluation of eq. (120) is 
cumbersome. It turns out that this is not necessary, if the time r 0 is chosen 
small enough. 

For elementary plane waves, the basic formula to justify this is the following : 



M fc r 0 < 1 


dy e ivy e i{k+x +k x+) , 


( 121 ) 


where k ± = (Mk/V 2) e ±y and = (tq/\/ 2) e. ±v . This formula can be checked 
by an explicit calculation of the integrals on both sides, which can be expressed 
in terms of Hankel functions for the right hand side and in terms of Gamma 
functions in the left hand side, and by doing the Taylor expansion to first order 
of the Hankel functions. The interpretation of this formula is the following: 

• the left hand side is the sum of the left- and right-moving partial waves, 
evaluated as if the time was r = 0 + (i.e. neglecting the evolution from r 
to r 0 ), 

• the right hand side contains the plane wave evolved to the time tq (i.e. 
the result of using the Green’s formula from r to tq). 
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In other words, eq. (121) shows that it is legitimate to neglect the time evolution 
of the spinors in the forward light-cone, provided that the time obeys M^tq <C 1. 
This exercise shows that if the time To obeys this condition, it is sufficient to 
add up the two partial waves on the light-cone, and to apply the transformation 
y —> v to their sum. Let us end this appendix by an important remark regarding 
the condition <C 1: it must be satisfied for all the transverse masses 

Mk = s/K + m 2 that can exist in the problem. In the lattice discretization of 
the Dirac equation, this implies that one must have To <C aj_ where a± is the 
transverse lattice spacing. 


D Conserved inner product 

D.l Definition and main properties 

Let us consider a locally space-like surface 16 E. For every point y £ E, we can 
define an orthogonal vector such that 

n^dy 11 = 0 for any displacement dy M around y £ E 
n° > 0 

= 1 . ( 122 ) 

Given two spinors i/’ and X, one can define the following inner product on E, 

= [ d 3 s v i>(y)ix{y)i ( 12 3) 

where d 3 S y is the 3-dimensional measure 17 on the surface E. One sees imme¬ 
diately that this inner product is Hermitean, 

(vix)* = (x|vO s • ( 124 ) 

The main property of the inner product defined in eq. (123) is that it is 
independent of the surface E if both ^ and x are solutions of the same Dirac 
equation 18 , 

(0)-m)il> = 0. (125) 

(This is true for any real valued background potential A^.) From now on, we 

can thus drop the subscript E in our notation for this inner product. Note also 
that this inner product is gauge invariant, since it involves the product of a 

16 This means that if is the local orthogonal vector to this surface, then > 0 . 

17 The measure on £ is defined in such a way that d 3 S y d(n ■ y) is the usual 4 -dimensional 
measure d^y. Therefore, there is some freedom in how we normalize the vector n^, provided 
we change accordingly d 3 S y in such a way that d 3 S y d(n • y) is left unchanged. 

18 Note that the quark spectrum involves the inner product between a spinor that has evolved 
over the background field and a free spinor. This inner product is therefore not conserved, 
reflecting the fact that the quark yield is time dependent and settles to a fixed value only in 
the limit r —y +oo. 
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spinor and the Hermitean conjugate of another spinor at the same space-tinre 
point. 

Let us give more explicit expressions of the inner product (123) for two 
important types of surface £. On a constant x° surface, it takes the following 
form, 

(dP |x) C onst *° = / d?, y i>Hy)x{y) ■ ( 12 6 ) 

On a surface of constant proper time, the integration measure is d 3 S y = Tdrjd 2 y± 
and the normal unit vector n ,J has the following components 


n = 


e +v 

V5 


n = 


o-fi 


V2 


n l = 0 . 


(127) 


Therefore, the inner product reads 


(^|x) T = t j dr)d 2 y± ^(r,y,y_ L ) e '' 7 ° t3 x(r,r/,y_ L ) . (128) 

D.2 Inner product at rr° = — oo 

This conserved inner product can be used as a consistency check for the various 
analytic formulas that we have obtained for the mode functions in the section 4. 
The first step is to evaluate the inner product on the surface y° = —oo , where 
the mode functions •0 fcsa are not yet modified by the background field. We get 

(27t) 3 2 uj k 8(k - k')5 ss '5 aa > , 

(27t) 3 2 u> k 8(k - k')5 ss 'S aa > , 

0 . (129) 

D.3 Inner product at r = 0 + in LC gauge and y basis 

Now, let us use the eqs. (68) and their counterparts for the left-moving partial 
wave in order to check that the spinors (in light-cone gauge, and in the y basis) 
evolved to the forward light-cone are consistent with eqs. (129). First of all, 
from the definition of eq. (123), we immediately see that 

(V’|x) r=0 + = J dy + d 2 y ± ip\y)V~ x(y) 

y-=0+,y+>0 

+V2 J dy~d 2 y±il)\y)V + x(y) ■ (130) 

y+=0+,y~>0 

In words, on the right branch of the light-cone we need only the V~ projection 
of the spinors, and their V + projection on the left branch of the light-cone. 


{^ksa life's' a') = 

i^ksal^k's'a') = 
(^ksal^k's'a 1 ) = 
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The other projections do not contribute to the inner product evaluated on the 
light-cone 19 . 

Adding the contributions of the two branches of the light-cone, we find the 
following expression for the inner product, 


(Vw|VWa') T=0 + = iV2(2ir) 2 5{k±-k , ± )5 l 


vi(k ± ,y)V + v s ’(k ± ,y') 


k+ - k'+ + ie 

v\(k ± ,y)T-v a >{k ± ,y'y 


k — k' + ie 


(131) 


where y, y' are the momentum rapidities corresponding to the 3-momenta fc, k'. 
Using Gordon’s identities, one sees that the imaginary part of the right hand 
side vanishes. Thanks to 

6(k + - k'+) = j^6(y - y f ) , S(k~ k'~) = j^S(y - y') , (132) 

where we denote M/~ = \Jk\ + to 2 , we arrive at 

0/’fcJV’fcw) r=o + = {2n) 3 S(y ~ y')6(k ± ~ k '^> 5aa 'JT k 

x vl(k±,y)(e~ v V + + e v V~)v s '{k ± ,y) .(133) 
The second line can be simplified by noticing that 

7 ° 7 3 =p+--p- , [V + ,V~]=0. (134) 

From these identities, we obtain easily 

e -yi°l 3 = e -yV+ e yV~ = e y v - + e -y v + ( 13 5 ) 

Using the fact that exp(— j/ 7 ° 7 3 /2) acts on spinors as a boost of rapidity — y in 
the z direction, we have 

vl(k±,y) [e v V~ +e~ v V + ] v s >(k±,y) = vl{k ± ,y)e~i J °' l3 e~^' 1 °" ,3 v s <{k ± ,y) 

= vl{k ± ,0)v s '(k ±1 0) 

= 2 M k 5 ss , . (136) 


Therefore, we obtain 20 

(V’fc S a|V’fc' s 'a') r=0 + = 2(2ir) 3 6(y-y')S(k ± -k' ± )6 ss ,6 aa , 

= 2w fc (27 T) 3 6(k-k')6 aa ,5aa> , (137) 

which is consistent with the conservation of the inner product. A similar verifi¬ 
cation can be performed for the positive energy mode functions ^ sa . 

19 Likewise, these projections do not contribute to the subsequent evolution of the spinors, 
because the relevant Green’s formula also contains a i/i. 

20 The second equality follows from S(y — y') = cu*. 5(k z — k' z ). 
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D.4 Inner product at r = 0 + in FS gauge and v basis 

We can perform the same check for the mode functions in the Fock-Schwinger 
gauge and v basis given by eq. (74). Firstly, let us apply the transformation 


'lpk±vsa — ^ dy e ^ ^k±ysa ( 13 $) 

to eqs. (129) in order to know what to expect for the inner product in the v 
basis. This trivial calculation tells us that we should obtain 


{i ) k J _vsa\ t l>k' ± v>8>a') = 2(27r) 4 (5(r' - is')6{k± - k' ± )5s S '5 aa ’ . (139) 

Let us now calculate this inner product directly from eq. (74). The inner 
product on a surface of constant r is given by eq. (128). It is convenient to 
absorb the factor r exp (— 777 °y 3 ) by defining new spinors that are the original 
ones boosted to the comoving frame at the rapidity 77 , times a factor yfr, 

i’(T,ri,yi.) = Vre ~^° l3 'ip(T,r],y±) . (140) 

In terms of these boosted spinors, the inner product reads simply : 

(V’|x) T = [ drjd 2 y± ^(r, 77 , y±)x{r, 77 , y±) . (141) 


When we insert two instances of the formula (74) in this equation, we see im¬ 
mediately that the crossed terms vanish because they contain (y + ) 2 = (y - ) 2 = 
0. Using the following identities, 


/ U 2(P±+ k ±) U 2(P±+k' ± ) 

ra + iv)ra-iv) = —^, 


= (27 r) 2 5(k ± -k' ± )S a 


(142) 


we first arrive at 

( ^k±vsa | i/'s'o') 


= (27r) 3 d(fcj_ — k!±)8(v — v')5 a 


Mfc COSh(7TI') 


xvt{k ±1 y = 0)\e™j y + + e ™y + y \v s ’(k±,y=0) . 

(143) 


The second line can then be simplified as follows : 


i’t(k±,y=0) e^y y + + e 7 ri, y + y 


= 2v\{k ± ,y = 0) 


e™V + + e 


- TTV ' D - 


'( k±,y = 0 ) = 
v s '(k±,y = 0) 


= 2vl{k±,y = 0) e™' f °' y3 v s >(k±,y = 0 ) 

= 2 nJ(fc_L, y = nv) v s >(k±,y=Tri') 

= 4 Mfc cosh( 7 ri/) S ss ' . 

Inserting this into eq. (143) gives the expected result (139). 


(144) 
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E Abelian case: electron spectrum in QED 

The initial value of the mode functions just above the forward light-cone, given 
in eq. (74), can also be used in the Abelian case. The analogous problem in 
QED would be that of the production of electrons in a high-energy collision of 
two large electrical charges Z\ and Z 2 . When they collide, the electromagnetic 
field of the two charges can produce electron-positron pairs. The spectrum of 
the produced electrons can be calculated by a formalism which is very similar 
to the Color Glass Condensate considered in this paper. In this description, the 
two colliding charges are replaced by the electrical currents they carry along 
their trajectories. These currents are the source of electromagnetic fields, that 
can be obtained by solving the Maxwell’s equations with sources : 


= Jl + J v 2 . (145) 

It is then this electromagnetic field that produces the charged fermions (in QED, 
this approach is known as the equivalent photon approximation.) 

Since QED is an Abelian gauge theory, it is much simpler than QCD in 
certain respects. Firstly, the Wilson lines U\ and U 2 are simply complex valued 
phases, that can be commuted at will. Secondly, each of the colliding charge 
produces transverse E and B which are localized in a shockwave transverse to its 
trajectory. In the Fock-Schwinger gauge, the corresponding gauge potential is a 
pure gauge in the half-space located after the shockwave. But since Maxwell’s 
equations are linear, the solution for the 2-charge problem is the sum of the 
solutions for individual nuclei, i.e. a sum of two pure gauge fields, which is itself 
a pure gauge. In QED, the fields are thus trivial inside the forward light-cone, 
unlike the QCD case. 

The formula (74), that gives the mode functions just above the forward 
light-cone, is therefore sufficient to obtain in closed form the amplitude. One 
can evaluate it by computing the inner product of eq. (10) at an infinitesimal 
time t —> 0 + , since the evolution at r > 0 is trivial. The only subtlety when 
doing so is that, since there is a non-zero pure gauge field in the forward light- 
cone, one should use a gauge rotated free positive energy spinor instead of 

= ul(x±_)Ul(x±_) u a (p) e~ ipx . (146) 

This free spinor must be first transformed as in eq. (70), 

^iw )+ W = VrU\ (x ± )ul(x ± ) J dy e wy e - ? W u<T (p) e _ip ' x . (147) 

In order to perform the integral over the momentum rapidity y, we need first 
to extract the y dependence hidden in the spinor u a (p ), 


u„(p±,y) = e ’* 3 


r(PJ_,0) = 


e*T + +e"ir 


u a (p _ L , 0) . 


(148) 
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Therefore, we have 


^S 2)+ (®) = yfru\{x ± )ul (x ± ) e ivr > 

x J dy e iuv e- iM * TCoMv) |e$P' +e ?'P ] ^(pj., 0 ) , (149) 

where M p = y^pj_ + to 2 denotes the transverse mass (in this equation, we have 
redefined the integration variable y—q—ty). The result of the integration over 
y can be expressed in terms of Hankel functions, thanks to 


dy e ay e 


Mv) = _ ine -i^ Hg\z) , 


where Ha\z) = J a (z)—iY a (z). In the limit r —> 0 + , we need only the beginning 
of the Taylor expansion of the Hankel function, 


H W(z) = i \(-Y a _-_ (-Y 

a y ; ^ 0 + sin(Tra) [\2/ T(1 - a) \2J 


a ' z-*o+ sin(7ra) [\2/ T(1 — a) V2/ T(l + a)J 

Note that when a has a non-zero real part, only one of the two terms domi¬ 
nates when 2 —» 0, depending on the sign of this real part. Therefore, in the 
combination s/zH^J ± (z), only one of the two terms survives when z —> 0 + , 


(z) = _ V2e^e^ 

tv 2 z->o+ cosh(7rzz)r(l — iv) 

Jze-^-^H^ Jz) = ^ e< * e ~ ¥ 

w +a z-> 0 + C 0 Sh( 7 TI/)r(| + iv) 


(r ■ <-> 


Therefore, 


$ { p±vv )+ { x ) t = Q+ ne U\{ Xl _)Ul{ Xl _) 


I 2 e lvri e ipj - <c± 
Mp C0Sh(7TI/) 




T (1 + iv) 


Similarly, the negative energy spinors evolved from the remote past up to r = 0 + 
read, 


^k,us( x ) = 


e ' 1 [ d 2 g_ l J_ 

fW k J ( 2 t r ) 2 M q 


M„T \ nfi , + ~ 

2M^J e 2 r (s -in)Ul(x ± )U 2 (q_L+k ± )'y-' 


e 2 r(i + in) Ul{x±)Ui(q± + k±) 7 “ 


x e lq± x± (q 1 7 * + m)v s (k ±1 y = 0 ) . 


38 



Using U + 7 + =V"f 
reads 


0 , the inner product between if) 


{U!U 2 ) + 

P± vo 


and 'I’hj.na 


(V& )+ |^J = - p) 


7re 2 
cosh(7r^) 


M p M k 


f d2q± - 1 t / n^ 

J (2tt)2 M q U ^ P±,0) 


U 2 {p±-q±)U 1 (q± + k ± )'y 

+ U 1 (p ± -q±)U 2 {q±+k ± )j + } ( 9 * 7 * + m) v s (k±, 0) . 

(155) 


In order to compare with existing results in the literature (e.g. eq. (52) in 
ref. [29]), one should perform the inverse transformation v —> y p ,p —> yk to 
return to momentum rapidity variables : 


« lC/ 2 )+ |V>J 


dvd d Jvyp P -ipyk 


(« 2)+ l^ 


k±fis 


)■ 


Thanks to the following formula, 


dv e ivz 
2n cosh(7r^) 


1 

27r cosh(|) 


(156) 


(157) 


it is easy to perform these integrals and one obtains 


(V^ lC/2)+ j K a )=iy/2 


d 2 q _l 

(2tt)2 


ut(p_L,y = 0) 


U 2 {p _ L - g_L)Ui(q_L + ki)e VP 2 Vk ^ 


+ 


M2 + 2 k~p+ 

M 

M2 + 2k + p~ 


Ui(p± - q±)U 2 (q± + fc_ L )e ! "' 2 yp 7 + 


■(m + g*7*)w s (fci,y = 0) . 

(158) 


The final step to compare with ref. [29] is to use the identities 


Vp—Vk + 

e 5 u1(p±,y = 0)j~(m + q l Y)v s {k ± ,y = 0) 

= ui(P±,yph~ (rn. + q l Y)v s (k ± ,y k ) 
y k -y P j. | 

e~^~ul(p±,y = 0 ) 7 +(?n + q l f)v s (k ± ,y = 0) 

= ul(p±,y P h + ('fn + q l Y)v s {k ±: y k ) , (159) 
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thanks to which we finally obtain 


(V£a lC/2)+ [^ fes ) = i ^f 4(P±,1/ P ) 

/ ^(px - g_L)^i(g± + fc-L)7~ 


V 


M2 + 2fc-p+ 


Ui(p±-q±)U 2 ( K q± + k ± )-y+\ {i 

- + 2fc+y . - (m + « 7 )«.(fex, ») • 


(160) 


This formula is identical to the eq. (52) in ref. [29]. Note that in order to 
recover this known result, it was crucial to gauge rotate the free spinor used in 
the projection, because the gauge field inside the forward light-cone is a non-zero 
pure gauge in the Fock-Schwinger gauge. 
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